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Abstract. We analyze a variational time discretization of geodesic calculus on finite- and certain classes of 
infinite-dimensional Riemannian manifolds. We investigate the fundamental properties of discrete geodesies, the 
associated discrete logarithm, discrete exponential maps, and discrete parallel transport, and we prove convergence 
to their continuous counterparts. The presented analysis is based on the direct methods in the calculus of variation, 
on F-convergence, and on weighted finite element error estimation. The convergence results of the discrete geodesic 
calculus are experimentally confirmed for a basic model on a two-dimensional Riemannian manifold. 

This provides a theoretical basis for applications in computer vision: In particular, discrete geodesies offer an 
effective tool for shape morphing and the computation of a distance between shapes, the discrete logarithm allows 
a linear representation of strongly nonlinear shape variability, the discrete exponential map provides a robust tool 
for shape extrapolation, and the discrete parallel transport can be used to transfer geometric details from one shape 
to another. The basic operations are exemplarily illustrated for two different shape spaces, the space of viscous 
volumetric objects and the space of discrete viscous shells. 
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1. Introduction. Riemannian geometry is a powerful theory which allows to transfer 
many important concepts (shortest connecting paths, the arithmetic mean, or the principal 
component analysis, to name but a few examples) from linear vector spaces onto nonlinear, 
curved spaces. It is based on the notion of a Riemannian metric, and it provides a set of 
basic tools useful in applications, among which we shall concentrate on geodesies, geodesic 
distance, exponential and logarithmic maps, and parallel transport. During the past decade, 
Riemannian concepts have for instance increasingly been applied in computer vision for the 
design and investigation of nonlinear and often infinite-dimensional shape spaces, where the 
Riemannian metric encodes the preferred shape variability. Applications include shape mor- 
phing and modeling [ ], computational anatomy in which the morphing path establishes 
correspondences between a patient and a template [ ], as well as shape statistics [ ]. 

Notwithstanding the conceptual power, operators such as the logarithm or the exponential 
map involve the solution of time-dependent nonlinear ordinary or partial differential equa- 
tions. In more complex spaces — as they appear in vision applications — they are typically 
difficult to compute unless the spaces possess very peculiar structures [31, 29]. Therefore, we 
here develop a discrete theory that can be seen as a natural time-discretization of the above- 
mentioned Riemannian calculus and which is comparatively simple to state and to implement. 
It is centered around the definition of a discrete geodesic as the minimizer of a time-discrete 
path energy and naturally extends from there to discrete analogs of logarithm, exponential 
map, and parallel transport. 

Complementarily to earlier work [ ], which deals with one particular shape space and 
focuses on the experimental study of the resulting discrete calculus, we here provide the 
theoretical justification of the general approach in terms of a rigorous convergence analysis 
for decreasing time step size under suitable assumptions on the manifold and the functional 
involved in the approximation. Furthermore, we discuss the range of applications and the 
method's current limitations. 

In what follows, let us give an overview on Riemannian calculus in the context of shape 
spaces in computer vision, which motivates our work. There is a rich diversity of Riemannian 
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shape spaces in the hterature. While some of them are finite-dimensional and represent shapes 
as polygonal curves or triangulated surfaces [17, 20], most deal with infinite-dimensional 
shapes such as planar curves which require Riemannian metrics with regularizing properties, 
for instance curvature-based or Sobolev-type metrics [ , ]. Interestingly, there is a strong 
link to mechanics. Based on the principle of least action, also complex mechanical systems 
like incompressible inviscid fluid flows can be rephrased in a Riemannian language, taking the 
kinetic energy as metric [ ]. Exploiting this idea, the space of sufficiently smooth domains 
(9 C can be assigned a Riemannian metric by identifying the tangent space at O with 
velocity fields v : O ^ and defining a (physical) metric on these. For instance, Dupuis et 
al. employ a metric g{v^v) = Lv -v dx for a. higher order elliptic operator L on some 
computational domain D C M.^ [ ], ensuring a diffeomorphism property of geodesic paths 
(cf. Section 3.3 below). In contrast, Fuchs et al. [ ] and Wirth et al. [ ] consider a second 
order elliptic operator based on viscous fluid flow, which is only defined on the shape domain 
O. The proposed discrete geodesic calculus might also be of use in optimal transport, where 
for two given probability densities po , Pi • ^ the Monge-Kantorovich functional 

J^\iIj{x)—x\'^Po{x) dx is minimized over all mass preserving mappings :D^D. Via a flow 
reformulation due to Benamou and Brenier [5], this fits nicely into the Riemannian context. 

There are only few nontrivial application-oriented Riemannian spaces in which geodesies 
can be computed in closed form (e.g. [ , ]), else the system of geodesic ODEs has to 
be solved via time stepping (e.g. [18, 4]). Geodesies can also be obtained variationally by 
minimizing the discretized path length [" ] or energy [i 1, 30], an approach generalized and 
simplified in this work. 

The idea of variational time discretization, which underlies our discrete geodesic cal- 
culus, has proved very appropriate also in other fields, in particular in the discretization of 
gradient flows for an energy with respect to a particular (Riemannian) metric (e.g. mean cur- 
vature flow [ ]). For Hamiltonian mechanical systems, variational time discretization is by 
now a classic field. The corresponding analog of our time-discrete path energy already oc- 
curred in the 70's as discrete action sum (cf. the brief historic account in [ ] ) and was 
analyzed as an independent, time-discrete system, whose discrete symplectic time steps ex- 
hibit the same structure as the discrete exponential proposed here. The consistency with the 
associated time-continuous mechanical systems was much later exploited numerically (cf. 
the general overview on variational, structure-preserving ODE integration [ ]) and analyzed 
for particular systems from the F-convergence ["^] and from the ODE-discretization per- 
spective [ ] under the name of variational integrators. These numerical schemes possess 
several advantages such as the (approximate) conservation of conserved physical quantities, 
the reproduction of the correct statistical properties, or the straightforward incorporation of 
additional system constraints, for which reason they have been applied in mechanical as well 
as computer vision contexts [12, 15]. More recently, time-adaptive versions have been de- 
veloped [ ], and the extremality of time-discrete action sums has also been applied as a 
constraint in optimal control [ ] . 

Let us note that in contrast to the above works, we here lay emphasis on the Riemannian 
tools, in particular on computing a discrete geodesic between two fixed given points, on 
deriving a consistent discrete logarithm from this geodesic, and on using discrete exponential 
and logarithmic map to compute a consistent discrete parallel transport. In Section 2 we 
briefly set forth the different components of our discrete geodesic calculus after which we 
show its applications to exemplary shape spaces in Section 3. The existence and convergence 
properties of the discrete geodesies and corresponding operators are proved in Sections 5 and 
6 after which we discuss how the theory applies on embedded manifolds and on certain shape 
spaces from Section?. 
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2. Discrete geodesic calculus on a Riemannian manifold. In this section we present 
the concept of the discrete geodesic calculus, which comprises the notions of discrete geo- 
desies, discrete logarithmic and exponential map, and discrete parallel transport. In prepara- 
tion for this we briefly recall the corresponding continuous calculus. Precise assumptions on 
the manifold and involved approximating functionals will be given in the subsequent sections. 

Let {A4^g) be a smooth, complete Riemannian manifold, where g denotes the metric, 
given as a family of positive definite quadratic forms gy : TyM x TyM M,for y e M, 
where TyM denotes the tangent space in y. Given a smooth path {y{t))te[o,i] on M, the 
length of this path is given by 

/:Mt))te[o,i]] = ^gy^t)m).m)^t. (2.1) 

Given two points yA and ys in the minimizer of C over all smooth paths (y(t))^^[o,i] 
with ?/(0) = yA and = is a geodesic. More generally, geodesies are defined as local 
minimizers of the path length. Let us define the distance dist(?/A, Vb) between ?/(0) = yA 
and y{l) = ys as the minimal path length. After reparameterization, a length-minimizing 
geodesic is also a minimizer of the path energy 

^Mt))te[o,i]]= f gyitMt).y{t))dt (2.2) 

and satisfies the constant speed property gy(t){y^ v) = ^'^[{y{t))te[o,i]]- The logarithm of a 
point ys with respect to a point yA is defined by log^^ ys = ^(0) if (^(t))tG[o,i] is the unique 
shortest path with constant speed connecting y{0) = yA and y{l) = ys- The exponential 
map exp^^ maps every tangent vector v G Ty^ M onto the endpoint of a geodesic starting at 
yA with initial speed v, i.e. in the above notation, exp^^ y{ff) = y{l) = yB • Furthermore, let 
us recall the parallel transport on a Riemannian manifold. Using the Levi-Civita connection 
V a tangential vector field v is parallel along a curve (?/(t))^^[o,i] on M if \/y(^t)v{t) = . If 
we sample a continuous path (?/(t))^^[o,i] at times = kr for k = 0, ... ^ K and ^ := 
denoting y^ = y(tk), we obtain the following estimates for length and energy 

K ^ K 

^[{y{t))te[o^]] > ^dist(?/fe-i,^/c) , ^[{y{t))te[o,i]] > -^dist^(?/fe-i,^/c), (2.3) 

k=l /c=l 

where equality holds for geodesic paths due to the constant speed property. Indeed, the first 
estimate is straightforward, and application of Jensen's inequality yields the second estimate, 

K K kr 

Vdist^(?/fe-i,y/c) < / gy(^t){y{t),y{t)) dt = r £[{y{t))te[o,i]] • 

k^i k=i ^(^-1)- 

Based on these preliminaries we are now in the position to introduce a discrete geodesic 
calculus. The starting point is a local approximation of the squared Riemannian distance 
dist^ by a functional W. In detail, we suppose that a smooth functional W : M x M ^ M. 
is given such that 

dist'(^, y) = W[y, y] + 0(dist^(^, y)) (2.4) 

for Sill y, y G M. We will see later that gy = ^W,22[^, y] implies (2.4) for smooth g and 
yV. Together with (2.3) this motivates the following definition of a discrete path energy 
and a discrete path length for discrete i^T-paths defined as {K + 1) -tuples (^o, • • • , ^k) with 
yk e MfoT k = 0, . . . ,K. 
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Definition 2.1 (Discrete length and energy). For a discrete K-path {yo^ - • - ^Vk) 
define the discrete length L and the discrete energy E by 

K K 

l[(?/o,...,^k)] = ^ \/y\^[yk-i,yk] , E[(yo,...,yK)] = K^y\}[yk-i,yk] • 

fe=l k=l 

Then a discrete geodesic (of order K) is defined as a minimizer ofE[{yo^...^yK)] for fixed 
end points yo^yx- 

To proceed with the definition of discrete logarithm and discrete exponential map we 
need to give a meaning to displacements y—y. Thus, in what follows let us assume that is a 
subset of some Banach space Y. This assumption is not restrictive because our constructions 
are purely local and can thus always be applied to charts of a more general manifold. Given 
a continuous constant speed geodesic (?/(t))^^[o,i] with y{0) = yA and y{l) = ys and a 
discrete geodesic (?/o, • • • , ^k) with y^ = yA and yx = ^b, we may view yi — yo as the 
discrete counterpart of ry{0) for r = Motivated by the fact that log^y^ {ys) = ^y{^) 
we hence give the following definition of a discrete logarithmic map. 

Definition 2.2 (Discrete logarithm). Suppose the discrete geodesic {yo^ . . . ^y^) is 
the unique minimizer of the discrete energy (2.1) with yo = yA cind yx = ys, then we define 
the discrete logarithm (-^LOG)^^ (?/b) = yi — yo • Note that is part of the symbol and 
not a factor 

In the special case K = 1 we obtain (-^LOG)^^ (ys) = ys — yA- As in the continuous 
case, the discrete logarithm can be considered as a linear representation (in the space Y of 
displacements on yA) of the nonlinear variation ys of yA- As we will verify later under 
suitable assumptions, for a sequence of successively refined discrete geodesies we obtain 
if(iLOG)^(y) ^ log,(y) fori^T ^ ^. 

In the continuous setting, the exponential map of a tangent vector v G Ty^M is defined 
via exp^^ v = ys, where ys = y{l) for the (unique) geodesic {y{t))te[o,i] with y{0) = yA 
and y{0) = v. Obviously, exp^^(-^i;) = y{^) holds for /c = 0, . . . , i^T. We now aim 
at approximating exp^^ (k-) via a discrete counterpart EXP^^ (the notation reflects the fact 
that on Lie-groups exp^(-) = exp(/c-)). The definition should be consistent with the discrete 
logarithm in the sense that EXP^^ (^i) = yk for a discrete geodesic (?/o, ^i, • • • , ^k) of order 
K > k with (i = yi — yo = (-^LOG)^^ (?/x)- Our definition will reflect the following 
recursive properties of the continuous exponential map for v e TyM sufficiently small, 

expy{lv) = {{ \ogy)~^ {v) , exp^(2v) = log^)"^ {v) , 

expy{kv) = exp^^p^((fc_2)^)(2v/e_i) for v^-i = log^xp^((;,_2)^) exp^((fe - l)v) . 

Replacing the tangent vector v by a displacement exp(/c-) by (EXP^), and ^ log by 
(^LOG) we obtain the following recursive definition of the discrete exponential map. 

Definition 2.3 (Discrete exponential map). For C ^ Y sufficiently small with y -\- ( e 
M define 

EXPI(C) = (iLOG);'(C) , EXP2(C) = (|LOG);'(C) , 

EXP^(C) = EXp2^p._.(^^(a-i) with Cfe-i = (iLOG)^xp._.(^)EXP^-i(C) . 

The smallness assumption on (" ensures that ( jLOG)^ and (^LOG)^ are invertible. It 
is straightforward to verify EXP?^ = (-^LOG) ^ on the image of (-^LOG) . The central 
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ingredient of this definition, the operator EXP^ requires the solution of a problem with a 
variational constraint. Indeed, to compute EXP^(C), one considers all discrete geodesic paths 
{y-, yi-,y2) of order 2, where for any chosen y2 the point yi must by definition be the unique 
minimizer of (2.1) so that we may write yi [y2]. Now, EXP^(C) is that point y2 for which we 
have ^ + C = ^1 [^2], which can be equivalently described via the equation 

?/ + C = argmin^^^_^ (>V[^,^i] + W[l/i,?/2]) . (2.5) 

The corresponding Euler-Lagrange equation, as a nonlinear equation for EXP^(C) = ^2, is 
given by W,2 [y, ^ + C] + W,i + C, ^2] = . 

Remark 2.4 (Discrete geodesies versus discrete exponential shooting). The varia- 
tional definition (2.5) of the operator EXP^ ensures that as long as a discrete geodesic 
connecting two end points yA, ys is unique, the discrete exponential map EXP^^ applied 
to (^LOG) {ys) for k = 0^ . . . retrieves exactly the same discrete geodesic as the 
minimization of the discrete energy E[7/^, ^1, . . . , yx-ij ys] over a// ^1, ... , yx-i^ i-^- yk = 
EXP^(;^LOG)^^ (ys) , where {yA = xo.yi, . . . , yx-i^ys = yK) is the discrete geodesic. 

Finally, we consider the discretization of parallel transport along a discrete (not neces- 
sarily geodesic) curve (?/o, • • • , yx)- In the continuous setting, given a curve (?/(t))^^[o,i] and 
a vector vq at ^(0), we define the parallel transport ^{y{t))t^^Q ij^o of vq along the curve as 
the vector v(l) resulting from the solution of the \/y(^t)^{t) = for t G [0, 1] and initial data 
v{0) = vq, where denotes the covariant derivative defined via the Levi-Civita connec- 
tion V. There is a well-known first-order approximation of parallel transport called Schild's 
ladder [ , ], which is based on the construction of a sequence of geodesic parallelograms. 
Given a curve (?/(t))^^[o,i] and a tangent vector v^-i G T^((fe-i)r)-^' the approximation 
Vk G Ty(^]^^^M of the parallel transported vector via a geodesic parallelogram can be ex- 
pressed as 

yl-l = exp^((/c-l)r)K-l) , yl = ^^'Pyl_, ' 

yl = exp^((/c-l)r) (2 log^((fe-l)r) (yl)) . = log^(/cr) (^D ' 

Here, y^ is the midpoint of the two diagonals of the geodesic parallogramm with vertices 
y{{k — l)r), yl_i, y^, and y{kr). This scheme can be easily transferred to discrete curves 
based on the discrete logarithm and the discrete exponential introduced before. 

Definition 2.5 (Discrete parallel transport). Let (^0, • • • , ^k) be a discrete curve in 
M with yk — yk-i sufficiently small for k = 1^ . . . ^K, then the discrete parallel transport of 
a displacement Co cit yo cilong (7/0, • • • , ^k) is defined for /c = 1, . . . , via the iteration 

yLi = yfc-i + Cfe-i , yl = vl-i + ((^log)^p_ Jyfc)) , 

yl = EXP2^_^ {yl - Vk-i) , Ck = vl - Vk , 

where is the transported displacement at y^. We denote Pyj^^...^yQCo = Ck- 

Above we have used (jLOG) (y) = y - y sls well as EXP^(C) = ^ + C- In the /c* 

1 y y 

step of the discrete parallel transport the Euler-Lagrange equations to determine y^ and C/c 
for given ^/c-i and discrete curve {yo^ • • - ^yx) are 

)A^Ayk-i + Ck-i.yl] + W,i [yl yk]=0, (2.6) 

w,2[^fe-i, + Wiiyl^yk + a] = . (2.7) 

If yV is symmetric, these conditions are the same as the Euler-Lagrange equations for back- 
wards parallel transport so that P^^,...,2yo ~ ^yo,---,yK^ which however is generally not true 
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if W is not symmetric. Note that for a continuous geodesic y{t) the velocity y{t) at a time t 
equals the initial velocity ^(0), parallel transported along the geodesic. In analogy, for a dis- 
crete geodesic {yo, . . . , yx) we have yk+i - yu = P^yfe,..., 2/0(^1 - ^0) for /c = 0, . . . , K - 1. 

Given discrete parallel transport, the Levi-Civita connection V ^r] for ^ TyM and a 
vector field r] in the tangent bundle TM can be also be approximated. 

Definition 2.6. For y e M, small enough ^ G Y, and vectors r]o attached to y and r]i 
attached to y -\- ^ (representing a discrete vector field), ^^{rjo^rji) = P~+^ y^i — ^0 defines 
a discrete connection. 

The inverse parallel transport satisfies the same Euler-Lagrange equations (2.6)-(2.7) as 
above (only this time y^ and (k-i are determined from (k)- In analogy to \/y(^t)y{t) = for a 
continuous geodesic a discrete geodesic (?/o, • • • , Vk) satisfies ^Ayk {^Vk^ ^Vk+i) = for 
^Vk = Vk+i — Ukr k = 0, . . . , iiT — 2. If W is symmetric, we can also express the discrete 
connection as V^{r]o, r]i) = Py^y^^rji - tjq. 

3. Examples for shape space manifolds. In this section we give some examples of 
Riemannian shape spaces used for shape analysis to motivate the discrete geodesic calcu- 
lus investigated in this paper. Currently the full analysis in this paper applies to the finite- 
dimensional models in Sections 3.1 and 3.2 and the infinite-dimensional model in Section 3.2 
as will be discussed in Section 7. 

3.1. The space of viscous-fluid objects. If we consider volumetric objects as shapes, 
i.e. open connected domains O C R^, we are in general not interested in point to point 
correspondences between two different objects. Hence, if we restrict to objects O which 
are homeomorphic to a reference object (9, then any shape y = O can be represented by 
an equivalence class of deformations [^] = | ^(O) = O}, and a path ((9(t))^^[o,i] of 
shapes can be represented by multiple families (^(t))^^[o,i] of deformations acting on the 
reference object O. A family {(j){t))te[o,i] is associated with a (Eulerian) velocity field v{t) = 
4>{t) o (j)~^{t), and shape variations are equivalence classes of such motion fields with v ^ v 
if V ' n = V • n on dO, where n is the outer normal on dO. One can define a physically 
motivated metric on shape variations as the minimal rate of dissipation in a Newtonian fluid 
occupying O when its free boundary moves according to the shape variation, 

go{VjV) = min / |^(tre[i;])^ + /itr(e[t;]^) dx , (3.1) 

{v \ v-n=v-n on dO} Jq 

where A, /i > 0, e[v] = \{Vv + Vv^). To approximate the squared distance between 
two objects O and O via a functional W, we solve a minimization problem over all possi- 
ble deformations between the two objects, i.e. >V[(9, (9] = min^^^^^^^^ dx for 
some particular energy density where A[il)] = Di/j^Dip is the associated Cauchy-Green 
strain tensor. To ensure the requested consistency (2.4) of the functional W and the metric 
we require W(l) = 0, DW{1) = 0, SiXid D^W{1){B, B) = ^{tiB)^ ^ i^tT{B^) for M 
B e R"^'"^ (cf. also Lemma 5.7 below). A suitable choice is W{A) = ^tiA + | det A - 
(/i + I) log Vdet A — ^ — ^ , which is rigid body motion invariant and isotropic so that 
W[QO + 6, Q(9 + 6] = W[0, O] for any two rotations and translations Q,Q e SO{d), 
b^b G R^. Hence, the resulting discrete path energy is also invariant with respect to rigid 
body motions independent of the time step size r (see [30, 26] for further details on this 
approach). Figures 3.1 and 3.2 show discrete geodesies between different shapes, where 
the colors indicate the local energy density W, which approximates the dissipation. The last 
three shapes in Figure 3.2 show extrapolations via the discrete exponential map. Figure 3.3 
depicts an example of discrete parallel transport. Both in the metric go and in the functional 
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Fig. 3.1. Discrete 8-geodesic between the digits 4 and 9 in the space of viscous fluid objects from Section 3.1. 
Colors indicate the local energy density W (from blue, low, to red, high), approximating viscous dissipation. 

3 3 3 77 7TV 



geodesic matching 



extrapolation 



Fig. 3.2. Discrete geodesic between the digits 3 and 7, extrapolated beyond 7, for the same model as in Figure 3.1. 



W the boundary of the objects is not sufficiently controlled to ensure existence of optimal 
continuous or discrete paths. As a remedy in what follows we restrict the admissible shapes. 
Alternatively, one might add a suitable regularization for dO (e.g. a perimeter regulariza- 
tion, [ ]) or combine the viscous-fluid approach for the object interior with the viscous shell 
approach (see Section 3.2) for its boundary. 

Following Fuchs et al. [a] the previous model can be reduced to a finite-dimensional 
shape manifold by describing the boundary of (9 as a B- spline curve with n control points 
qi, . . . ^Qn. Hence, each shape is represented by the vector y = (gi, • • • , ^n)- This time 
we pose the shortest path problem between two B-spline curves under the assumption of 
a point correspondence between these curves. The evaluation of the metric thus requires 
to solve the associated minimization problem (3.1) for given velocity of the control points 
^ = y = {Qi^ • ' ' ^Qn) (and thereby of all points on the curves). Instead of the nonlinear 
function W above we might consider for two sets of control points y and ^ as in [ ] the 
quadratic functional yV[y^y] = Jo^y^ |(tre['^[?/, ^]])^ + /itr(e[?/^[?/, ^]]^) dx where 0[y] is 
the object domain described by y and ^/^[t/, y] is the energy-minimizing deformation under the 
Dirichlet boundary conditions given by the point to point correspondence between the two 
B- splines. 



3.2. The space of viscous shells. Here we deal with shapes S which are shells of some 
small thickness S > 0. Deforming these thin material layers, the associated viscous friction 
is predominantly caused by two mechanisms, namely friction due to tangential (trans versally 
uniform) in-plane deformation and friction due to (transversally non-uniform) deformation 
caused by bending. Assuming a point to point correspondence between two shells S and 
iS, similarly to before we can define their squared geodesic distance as the minimum energy 
dissipated during a distortion of S into S. It can be approximated by a functional measuring 
the strain energy of a deformation (j) : S ^ with (j){S) = S. The tangential strain is 
described by the tangential Cauchy-Green strain tensor A [(/)] , which is defined pointwise as an 
endomorphism on the tangent space of 5 wisi I {A[(j)]v , w) = I{d(j){v)^ d(j){w)), where / and 
/ are the first fundamental forms on S and 5, respectively. As the corresponding energy we 
apply the same functional as in the volumetric cases above, this time acting on the tangential 
Cauchy-Green strain tensor, i.e. >Vtgi[y, y] = S W{A[(j)]) da for (j) = y o y~^ , where y, y 
are parametrizations of S and S, respectively, over a reference shell S. The dissipated energy 
due to bending can be approximated as follows. Let S, S denote the shape operators on S 
and S, respectively. Furthermore, define S[(j)] as the pullback of the shape operator S onto S 
given by w) = I{Sd(j){v)^d(j){w)). Now the bending energy of the deformation (j) 



8 



MARTIN RUMPF AND BENEDIKT WIRTH 



discrete geodesic — , 




Fig. 3.3. The digit 4 is sheared or bulged, and the shearing and bulging variations are transported along the 
discrete geodesic to the digit 9 via discrete parallel transport based on the model from Section 3.1. 



Fig. 3.4. Discrete 8-geodesic between two hand poses, based on the model of discrete viscous shells from 
Section 3.2. The colors indicate dissipation due to surface stretching (top) and bending (bottom). 

can be described via the difference Q[(j)] = S[(j)] — S, known as the relative shape operator, 
and it is expected to scale cubically in 5 (cf. the rigorous Hmit analysis i n [ ]). F or the sake 
of simplicity we choose Wbnd[^7 v] = Is \Q[^]\'^ ^ where \Q\ = y^tiQ^Q denotes the 
Frobenius norm. Combining both energy contributions we finally obtain 



with (j) = y o . To apply this approach in geometry processing we consider triangu- 
lar meshes as discrete approximations of shells, where each shell is represented by a vec- 
tor = (xi)i=i,...^n ^ E^"^^ of n nodes Xi G R^. The gradient of a deformation 
(j) = y o y~'^ is piecewise constant on each shell triangle so that is also piecewise 
constant. To define a discrete shape operator on a triangle T we proceed as follows. For 
every edge e = Ti fl T2 between two triangles, we obtain an edge normal Ne as the nor- 
malized sum of the unit normals on Ti and T2. Then, for a triangle T with cyclically 
ordered edges eo, 61,62 and corresponding edge normals No^Ni,N2 we define the dis- 
crete shape operator as the piecewise constant endomorphism S on each triangle plane with 
I{Sek, 6/e) = 2(7V(/,+i) mod 3 - ^(fc-i) mod 3)e/c. From these discrete shape operators we 
then derive the piecewise constant relative shape operator as above. Finally, we obtain a dis- 
crete counterpart of the energy (3.2) as the basis of a time-discrete geodesic on the space of 
triangulated shells. For details on this approach we refer to [ ]. Let us remark that this 
approach can be combined with the shape space model for viscous-fluid objects, where we 
simultaneously measure viscous dissipation on the object domain and its boundary contour. 
Figure 3.4 shows a discrete geodesic between two shells, color-coded with the local energy 
density in (3.2) due to both tangential and bending deformation. 






(3.2) 
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If we restrict ourselves to one-dimensional shells in and parametrize them over the 
interval [0, 1] we obtain 

Jo \\y,s\ J Jo 

for the tangential distortion and the bending energy, where K:[y] = (^^) ' ^\y T 
curvature of a curve y, subscript ,s denotes the derivative, and D^^ counterclockwise rotation 
by |, i.e. n[y] = ^| is the normal on the curve y. The underlying Riemannian metric gy 
on the space of shell curves for which ^ W,22 [y^ y] = 9y is given by 

gy{v,v) = f 25^W"{l)\vf\^ + 5^ {dyn[y]{v)f As , 

^0 \yM 

where v^f = ' denotes the tangential component of v^s and 

\y,s\ ' \y,s\ \\y,s\ \ \y,s\ \y,s\J 

Linearization of the bending energy and the choice = ^"^"^^ lead to a simplified model 

W[y,y] = ^'^(i-|^) \y,s\ + S'{y,ss-y,ss?\y,s\ds, (3.3) 

which corresponds via gy = ^ W,22 y] to the simplified metric 

pi \ytg\\2 

gy{v,v)= 2S'-f^^S^\v,ssf\yAds (3.4) 

Jo \y,s\ 

to which our convergence theory will apply as discussed in Section 7. 

3.3. The flow of diffeomorphism approach. Here we consider a domain 1] c for 
(i = 2,3, and families of diffeomorphisms (0(t))^^[o,i] : ^ H^. In medical applications 
[ ] the diffeomorphisms represent deformations of anatomic reference structures. Thus, each 
diffeomorphism (j) : Q ^ M.^ represents a particular anatomic configuration or shape y = (j) 
in the space M of all such diffeomorphisms. We associate the energy ^ [((/>(^))tG[o,i]] = 
Lv{t) ■v{t) dxdt with a path (^(t))^^[o,i] in this shape space, where v{t) = <j){t) o 
(t){t)~^ represents the Eulerian velocity of the underlying flow and L is a higher order elliptic 
operator. Under suitable assumptions on L one can prove that paths of finite energy are indeed 
one-parameter families of diffeomorphisms as studied by Dupuis, Grenander and Miller [ ]. 
Physically, the metric g^j^^t) ^(^)) = Jq Lv{t) ■ v{t) dx describes the viscous dissipation 
in a multipolar fluid model as investigated by Necas and Silhavy [ ] . A specific model for 
the viscous dissipation is given by (j)) = |(tre['u])^ + /itr(e['u]^) + j\D'^v\'^ dx , 
X^jJ^jj > 0,m > 2. The first two terms of the integrand again represent the usual dissipation 
density in a Newtonian fluid, whereas the third term represents a higher order measure for 
friction. Obviously, this model is invariant with respect to rigid body motions. 

A suitable choice for the functional W approximating the associated squared distance is 
given by 4>] = Jq W{A[ijj])) + 'y\D'^ijj\'^ dx for = o (/)~^, where W is the energy 
density from Section 3.1. Once again, the resulting discrete path energy would be invariant 
with respect to rigid body motions. 
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4. Discrete geodesic calculus on embedded finite-dimensional manifolds. After hav- 
ing introduced a number of relatively complex shape manifolds as a motivation for its ap- 
plication, let us briefly demonstrate the discrete geodesic calculus for the much simpler case 
of an (m — 1) -dimensional manifold A4 embedded in R^. We consider the simple energy 
y^[y^y] = \y~y\'^ which reflects the stored elastic energy in a spring spanned between points 
y and y through the ambient space of M. in R^. Thus, the discrete path energy is given by 

e[(^o,...,^k)] = ^Ef=i l^/c -^fc-ip. 

Now, we define the Lagrangian F(X, A) = E[(?/o, • • • — A • D{X) , where X = 

(^1, • • • , yK-i) are the actual point positions, A = (Ai, . . . , Xr-i) the Lagrange multi- 
pliers, and D{X) = {dM{yk))k=i,...,k-i with dM{') the (local) signed distance function 
of M. Hence gradj^F(X, A) = 0, grad^F(X, A) = are the necessary conditions for 
(t/o, . . . , ^k) to be a discrete geodesic with end points y^ = yA and yx = ys- The resulting 
system of (m+ 1) (i^T — 1) degrees of freedom can be solved for instance by Newton's method. 
In particular, (gradj5^F(X, A))k = 2K {2yk - yu-i - yu+i) - Xk^dMiyk) = results in 
finding points yk (k = 1, . . . , K — 1) on M such that 

2yk - yu-i - yu+i ^ Ty^M , (4.1) 

while grad^F(X, A) = D{X) = ensures that the points y^ are located on M. Equation 
(4.1) is a discrete counterpart of the continuous geodesic equation y{t) ± Tyt^^^M. for an 
arclength parametrized geodesic curve {y{t))teR with '^Vk-yk-i-yk+i ^ _y(^]^^y 

From these considerations, we derive for the discrete logarithm that C = (|LOG)^y is 
described by the condition 

To achieve a simple geometric condition for the discrete exponential map EXP^^ we consider 
a discrete geodesic (?/, ^i, ^2) with yi = y ^ ( and ^2 = ^1 + where the displacement r] is 
the actual degree of freedom. The necessary condition for (?/, ^1,^2) to be a discrete geodesic 
and thus for y2 = EXP^^ is 

C-r]± Ty^^M , 

which results in a one-dimensional search problem in the space fl (?/ + span{C, n^y+^j), 
where Uyj^c, denotes a normal on Tyj^c^M. 

5. Properties of the discrete path energy. In the remainder of the paper we aim to 
examine the convergence properties of the discrete geodesic calculus as the discrete time step 
size T = ^ tends to 0. Figure 5.1 shows experimental evidence for the convergence of the 
discrete geodesic, exponential map, logarithmic map, and parallel transport. In this section 
we study some properties of the discrete and the continuous path energy on a specific class 
of manifolds, and we examine their relation via the concept of F-convergence. In particular 
we will show that sequences of successively refined discrete geodesic paths converge to a 
continuous geodesic path. 

We consider the following functional analytic set up. Let Y be a Banach space with norm 
II • II Y on which a symmetric bilinear form (?/, w) ^ gy{Vj w) for y^ w e Y induces a 
Riemannian structure with ^:YxYxY^IRU{oo}. Furthermore, let V be a compactly 
embedded, separable, reflexive subspace of Y with norm || • ||v. We shall suppose that g is 
uniformly bounded and V-coercive in the sense 



c*\\vf^<9y{v,v)<C*\\v\\l 



(5.1) 
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Fig. 5.1. Top: Discrete geodesic path yr = (yo, • • • , yx) ^ith yo = yA = (§5 0) cindyi = ys = ("I? 2) 
on the stereographic projection of the sphere, X(^LOG) (EXP'^) (v/K) for v = log (ys), (ind 

parallel transport ofw = ( — | , 0) from yA to yB at different temporal resolutions (K = 2^ for k = 1, . . . , lOj. 
Bottom: The corresponding errors versus 1 /K. 



for constants c* > and C* < oo independent of y. In particular, gy{v, = oo if v V. 
Also, g shall be continuous in y with 

\9y{v,v)-gy{v,v)\</3{\\y-y\W)\\vf^ (5.2) 

for a strictly increasing, continuous function /3 with /3{0) = 0. Finally, given G Y we 
consider the manifold A4 := + V, on which ^. (•, •) induces a Riemannian structure. 

Remark 5.1. As a motivation for this setup, consider the infinite -dimensional shape 
space of ID shells discussed in Section 3.2 ( the application to finite -dimensional shape spaces 
is explained in the next remark). Shells y can be considered as differentiable curves in the 
space Y = C^([0, 1]; R^), and the space V = VK^'^((0, 1); R^) encodes the regularity of 
admissible variations vofa shell y for which the metric gy{v,v) is finite. The shell manifold 
M is locally spanned by all those shapes which can be reached from a particular ^ Y 
via (short time) integration ofy(t) = v(t) with v(t) G iy^'^((0, 1); R^). 

Remark 5.2. The assumptions particularly cover the local theory of smooth k-dimen- 
sional shape manifolds (k < oo), which by Nash's theorem can be isometrically embedded 
in the Euclidian space R^ for m sufficiently large. In this case one chooses Y = V = R^ 
to represent a chart of the manifold and gy{v,v) = DX{y)v • DX{y)v, where X is the 
associated parametrization and • is the scalar product in R^. 

For yA^ yB ^ A1, the next theorem states the existence of a connecting path with least 
energy. The key point is the weak lower semi-continuity of the continuous path energy (2.2) 
using the compact embedding of V into Y. 

Theorem 5.3 (Existence of continuous geodesies). Under the above assumptions there 
exists a continuous geodesic path (?/(t))^^[o,i] defined as a minimizer of the path energy £ 
among all paths (^(t))tG[o,i] "^i^h y{0) = yA andy{l) = ys- 

Proof Let {y^ {t))t^^Q^i\ be a minimizing sequence of paths. We can obviously assume 

that the energy on this sequence is bounded by f = gyj,+t{yB-yA) (Vb -Va^Vb- Va) dt. 
Hence, ||l2((o,i);V) is bounded for = y^ . Because i^^((0, 1); V) is separable and re- 
flexive, there is a subsequence — with a slight misuse of notation indexed as before — which 
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converges weakly to some v G L^((0, 1); V). Likewise, the curves {t) = v^{s) ds con- 
verge weakly in iy^'^((0, 1); i/a + V) to y{t) = v{s) ds. Due to the bounded embedding 
of V^i'2((0, 1); V) into C^'^ ([0, 1]; V) and the compact embedding of C^'^ ([0, 1]; V) into 
C^([0, 1]; Y) by Arzela-Ascoli, upon selecting a further subsequence, converges strongly 
in C^([0, 1]; Y). Now choose j large enough that E[y^] < E-\- S and sup^^[o,i] PiWv^i'^) ~ 
^(^)||y) < for a prescribed small S > and S = infS. Following the usual paradigm 
of the direct methods in the calculus of variations we can estimate the energy using Mazur's 
lemma, the convexity of v gy{v,v), the continuity and coercivity of g, and Fatou's lemma, 

£Mt))te[o,i]] = I gyitMt),v{t))dt 
Jo 

< / lim Y.}^^gy^,^{v\t)y{t))dt 
Jo ^^-t^ 

Here, (A^j)j=i,...,oo is a sequence in IN and (A;^)^"^^'"'^ is a sequence of convex combination 

coefficients with Yjf=j K = 1' K — 0' ^i^j -^i^* ~^ ^ strongly in i^^((0, 1); V) for 
j ^ oo. Because S is arbitrary, we obtain £[{y{t))te^Q^i]] < S, which proves the claim. □ 
For yAj ys ^ Y, let us introduce the Riemannian distance 

dist(yA, ys) = . / min £[{y{t))te[OA]] • (5.3) 
Y yio)=yA,yii)=yB 

It is an easy exercise to verify the axioms of a metric and that the induced topology is equiva- 
lent to the V-topology, Vc^\\yB-yA\W < dist(yA,te) < VC^\\yB-yA\W- Furthermore, 
a simple reparameterization argument shows gy(^t^{y{t)^y{t)) = const, along a geodesic. 

For the discrete path energy, we would like to show an analogous existence result, as 
well as properties related to (5.3) and the above-mentioned constant speed parameterization. 
For this purpose we consider a lower semi-continuous W:YxY— ^IRLJ{oo} satisfying 

W[y, y] = dist2(y, y) + 0(dist3(?/, y)) (5.4) 

with uniform constants, and we assume coercivity of W in the sense 

W[y,y] >7(dist(^,^)) (5.5) 

for a strictly increasing, continuous function 7 with 7(0) = 0. 

Theorem 5.4 (Existence of discrete geodesies). Given yA-, ys ^ M.^ there is a discrete 
geodesic path {yo^ • • • ^yx) which minimizes the discrete energy E over all discrete paths 
(^0, . . . , yx) with yo = yA and yx = ys- 

Proof. Let ((?/q, . . . , ^;^))j=i,...,oo be a minimizing sequence. We can obviously as- 
sume the yl, to lie on = + V and that the energy on this sequence is bounded by 
E = KyV[yA-,yB]^K{K—l)W[yB-,yB\' Because V compactly embeds inY, there is a sub- 
sequence, still denoted {y^^ - - - ^y^)^ which converges in Y to a discrete path (^0, • • • , ^k). 
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By the lower semi-continuity of W in both arguments, we finally obtain E[(?/o, • • • , ^k)] < 
in%o,...,yx)E[(^o,...,^K)]. □ 

Theorem 5.5 (Bounds on discrete path energy). The discrete path energy satisfies 

min E[(?/o,...,yK)] =dist^(yA,yB)(l + 0(dist(yA,^B)/i^)) 

(?/Ov,?/k) 

yo=yA,yK=yB 

with uniform constants. 

Proof. Denote the minimum energy value by E. Choosing points ^i, . . . , yK-i such that 
dist(?//c_i, = dist(?/o, yK)/K for k = 1, . . . , we directly see 

E < i^Ef^i rnVk-uVk] = dist\yo,yK){l + 0{dist{yo,yK)/K)) . 

On the other hand, letting ak = dist{yk-i, yk) for the minimizing (?/o, • • • , ^k), we have 
E>K Ylk=i ~ ^^k ~' • • • 5 ^k) for some C > 0. From the coercivity of W we 

knowE > KEf=i7K) and thus afc < -f-\dist\yo,yK){^^0{^^^%^))). Thus, for 
K large enough, ak < 1/3C for all k. Minimization of F under the constraints EfeLi > 
dist{yo^yK) and ak < 1/3C, /c = 1, . . . yields ai = . . . = = dist(?/o, y^j/i^ and 
thus E > F(ai, . . . , ax) = dist^(yo, - C'dist(?/o, yK)/K). □ 

Theorem 5.6 (Equidistribution of points along discrete geodesies). Discrete geodesies 
( minimize rs from Theorem 5.4) satisfy dist{yk-i^yk) ^ Cdist{yA^ yB)/Kfork = 1, . . . , iiT . 

Proof. For given K, let Jk ^ {1, . . . , i^} be such that •= dist(?/j^_i , ) is largest. 
Furthermore abbreviate d := dist(?/o, yx) and := dist(?/o, have 

d2(l + 0(d/if)) = E[(yo,...,yif)] 

> 7]^E[(yo, . . . , Vj^-i)] + K^l{l - 0(7if)) + K^niVm Vk)] 
>if7l:(l-0(7K)) + j^dist2(yo,%K-i)(l-0(7K)) + 7^dist2(y,-,,y^)(l-0(7K)) 

> (1 - 0{^k)){K^I + + K^{l-aK- ^-f?]) . 

This is minimized by aK = — ^ and yields 

d'{l + 0{d/K)) > (1 - 0{^k)){K^I + - ^)') . 

For 7k = ^ with ^ oo for K ^ oo and ujk < Ca/K by the trivial estimate < 

3 2 

one obtains + ^ > + ^^l^i^ ~ )^) ^^^^^ ^^^^^^ ^ contradiction, 

so dist(?/^-^_i,?/^-^) < Cd/K. □ 

For the subsequent estimates it is convenient to relate the function W to the metric g; we 
use the following. 

Lemma 5.7 (Consistency conditions). Ifv ^ y-\-v] is twice Gdteaux-differentiahle 
on \ for y gY, then W[y^ y] = dist^(?/, y) + 0(dist^(?/, y))for y close to y implies 

W[y,y] =0, W,2[y,y](^) =0, W^22[y.y]{v,w) = 2gy{v,w) 

for any v^w G V. Furthermore, W,i y] iy) = and 

W 11 [y, ^] ^) ->V,i2 ^] {v, w) = - W,2i y] {v, w) >V,22 [y, y] {v, w) . 

If V yV[y^y -\-v] is even three times Frechet-differentiable, the implication becomes an 
equivalence. 
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Proof. It is readily shown that dist^(?/,^) = gy{y-y,y-y)^0{j3{\\y-y\\Y))\\y-y\\\i- 
Let W[^, y] = dist^(?/, y) + 0(dist^(?/, y)) for y close to then obviously W[y^ y] = ^ 

and W,2[^, ^]('^) = + = 0- Now assume W,22[^, ^]('^, 7^ 2^2y('^5 f^^' ^^^^ 

V G V (note that due to the bilinearity it is sufficient to show equality of the quadratic forms). 
Without loss of generality, let W,22[^, ^) > agy{v^v) for some a > 2. This implies 
yV[y^ y -\- tv] > ^ (f + 1) 9y{tv^tv) for t small enough, which (using dist^(?/, y -\- tv) = 
gy{tv^ tv) + o(t^)) is strictly greater than (^ + ^)dist^(?/, y -\-tv) -\- o(dist^(y, y + t'u)) and 
thus a contradiction. 

The above applied to the first argument of W instead of the second implies W,i [y, ^] =0 
and yV^ii[y^y]{v,w) = 2gy{v,w). Finally, for any curve {y{t))teii in Y we can differen- 
tiate = W^i[y{t),y{t)] and = W^2[y{t),y{t)] with respect to t, yielding W,i2[^,^] = 
-W^ii[y,y] and W,2i[^,^] = -W,22[^,^]. 

If V W[y ^ y-\-v] is three times Frechet-differentiable and W,22 ^] ^) = ^^gyi^^^)^ 
then by Taylor's theorem foYy = y-\- tv, W[y^y] = gy{tv^tv) + 0{t^) = dist^(?/,^) + 



Theorem 5.8 (Uniqueness of discrete geodesies). Ifv ^ yV[y^y-\-v] is twice Gdteaux- 
dijferentiable on W for y G A4 and if\\yB — ^a||v is sufficiently small there exists a unique 
discrete geodesic {y^^ . . . ^yx) y^ith yo = yA crnd yK — yB- 

Proof Without any restriction we assume ^ V. Otherwise, we have to work with 
offsets instead of points For Xo,K := (^o, ^k) and := {yi, . . . ,yK-i) 

we define the function F : M^+i ^ (V^^-^ (Xo,k, ^ (W,2[yfe-i, ^fe] + 

W,i [yki yk+i])k=i,...,K-i , where Y' denotes the dual space of V. Then the block tridiagonal 
operator A = I^Xi ...,k-i^(-^o,k, Xi,...,k-i) is sl {K — 1) x {K — 1) is given by 



where g = gy^ . Hence, the inverse of A can be computed by Gaussian elimination, which 
ensures that A = • • • , ^a) is invertible with bounded inverse, and thus the 

claim follows by the implicit function theorem. □ 

In what follows, we aim to prove convergence of discrete geodesies against continu- 
ous ones. To this end we identify any discrete path {yo, . . . ^yx) on M. with its piecewise 
geodesic interpolation y^ : [0, 1] M, i.e. every segment ^^|[/cr,(fe+i)r] shall be the short- 
est continuous connecting geodesic between y^ and yt+i- Now, we define an energy 
on continuous paths via ^^[{y{t))te[o,i]] •= E[(yo, • • • , ^k)] if y{t) = y^ {t) for some 
yo,...,yK ^ M with yo = yA, yK = ys, and ^^[{y{t))te[Q,i]] = ^ else. Based on these 
notational preliminaries we obtain the following convergence result. 

Theorem 5.9 (F-convergence of the discrete energy). In the i^^((0, 1);Y) -topology 
the V -limit of for K ^ 00 is £ . 

Proof. To verify F-convergence we have to establish the two defining properties, the 
limsup- and the liminf-inequality. 

[limsup -inequality] For an arbitrary : [0, 1] ^ with y G i^^((0, 1); Y) we have to show 
that there exists a sequence of curves y^ : [0, 1] ^ with y^ ^ y in i^^((0, 1); Y) and 



0{disi\y,y)).U 



Akk = W,22[yfe-l,^fe] +>V,li[y/c,^/c+l] , 
A(/c-l) = W,2l[yfe-l,^fe] , Ak{k+1) = W,12[^fe,^fe+l] . 



Now, using Lemma 5.7 we obtain for the trivial geodesic {yAj - - - ^yA) that 
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limsup^^^ [{y^ {t))te[o,i]] ^ ^[(^(^))^g[o,i]] • Without any restriction we can assume 
^[(^(^))tG[o,i]] < CO. Applying the Cauchy-Schwarz inequahty we get 

dist{y{t), y{s)) = j^\J 9y{r) {y{r) , y{r)) dr < yJ\s-t\J' Qy^r) {y{r) , y{r)) dr , (5.6) 

which immediately implies iy^'^((0,l);A^) C C^'^([0,1];A^) and thus Holder continuity 
of paths with finite path energy. Now, let y^ denote the piecewise geodesic interpolation of 
(^(^), • • • , y{§))- For any K we have 

niymeio,!]] > i^Ef=idist^(y(^),^(#)) 

>i^Ef=iW[^(^),^(#)]-Gi^Ef=idist^(^(^),^(#))) 

> E^[(^^(t)),,[o,i]](l - CK-iy^Sm))^]) , 

where we have used (5.6). Letting K ^ oo yields the desired limsup-inequality. 
[liminf- inequality] We have to show that for any sequence of curves {y^)K>i with : 
[0, 1]^M and y^ ^y'm L\{0, 1); Y) the inequality liminf^^oo E^[(^^(^))^g[o,ii] > 
^[(^(^))tG[o,i]] holds. Without any restriction we assume that [{y^ {t))te[o,i]] < E < 
oo uniformly. Thus, yV[y^ {^^)^y^ {^)] uniformly as i^T ^ oo so that due to the 
coercivity dist(^^(^^), ?/^(-^)) uniformly as well. Next, we estimate 

hminf E^[(y^(t)),,[o,i]] = liminf /fEf=iW[y^(^),y^(|)] 

K^oo K^oo 

> liminf ifEf=idist^(y^(^),t/^(#)) (1 - Cdis%^(^),y^(A ))) 

= liminfE[(y^(t))^e[o.i]], 

which also shows that y^ is uniformly bounded in iy^'^((0, 1); A^). Due to the reflexivity 
of V, a subsequence (for simplicity again denoted by {y^)K>i) weakly converges against y 
in VK^'^((0, 1);M). Due to the sequential weak lower semi-continuity of the energy E (cf. 
Theorem 5.3) we finally obtain the requested estimate 

liminf E^[(y^(t)),e[o,i]] > liminf £:[(y^(t)),g[o,i]] > £Mt))te[o,i]] , 

which concludes the proof. □ 

Corollary 5.10 (Convergence of discrete geodesies). Minimizers of the discrete path 
energies E^ converge against minimizers of the continuous path energy £ in C^([0, 1]; Y). 

Proof This is a simple implication of the F-convergence and the following equi-mild co- 
ercivity of the discrete energies: Theorem 5.5 shows that the minima of the discrete energies 
are uniformly bounded. However, E^ is coercive with respect to yA + I^^'^((0, 1); V) so 
that the minima for all K are achieved in some bounded ball of yA + W^^'^ ((0, 1); V), which 
is compact in L^((0, 1); Y). Together with the F-convergence of E^, this implies that any 
L^((0, 1); Y) -limit point of a sequence of minimizers y^ for E^ is a minimizer for £. Fur- 
thermore, from the above, all converging subsequences are bounded in yA + W^^'^((0, 1); V) 
and thus bounded in + C^' 2 ( [0, 1] ; V) and by Arzela-Ascoli precompact in ( [0, 1] ; Y) 
so that the convergence is not only in i^^((0, 1); Y), but even in C^([0, 1]; Y). □ The above 
convergence obviously does not only hold for the piecewise geodesic interpolation y^ of dis- 
crete geodesies (yo, • • • , ^k), but also for the piecewise linear interpolation which we shall 
call yr (where r = 1/K stands for the discrete time step). 
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For Stronger convergence estimates and for the convergence of discrete logarithm, expo- 
nential map, and parallel transport, we require the following smoothness hypotheses: 
(HI) The metric g is C^(Y; V (g) V')-smooth. 

(H2) The energy W is + V) x {y^ + V); IR)-smooth with bounded derivatives. 

The following theorem now states that the convergence in C^([0, 1]; Y) ensured by the 
above F-convergence result is actually much stronger with -converging velocities. 

Theorem 5.11 (Path convergence in W^'2((0, 1);M)). Under the hypotheses (HI) 
and (H2), if the interpolated geodesies yr{-) converge for K ^ oo to the continuous one y{-) 
in i^^((0, 1); Y), then this convergence is even in W^'^((0, 1); Ai). 

Proof Under the smoothness and coercivity assumptions on g, a continuous geodesic 
path (?/(t))^^[o,i] in as a minimizer of the continuous energy gy(^t^{y{t)^y{t)) dt for 
given fixed end points y{0) and y(l) is smooth and fulfills the Euler-Lagrange equation 



for all G VFo'^((0,l); V) := {V^ G VFi'2((0, 1); V) | V^(0) = V^(l) = 0}. Furthermore, the 
Euler-Lagrange equation for a discrete geodesic path is given by 



for all {il^k)k=o,...,K C V with ?/^(0) = V^(i^) = 0. Applying the Taylor expansion 

W[yk-i,yk] = yv[yk-i^yk-i] + >V,2[^/c-i,^/c-i](^/c - Vk-i) 

+ / {i-s)W^22[yk-i^yk-i^s{yk-yk-i)]{yk-yk-i^yk-yk-i)ds, (5.9) 




(5.7) 



KELi i^Avk-i^ykK^k-i) ^W2[yk-i.yk]{^k)) = o 



(5.8) 



whose first two terms on the right-hand side vanish, one can rewrite (5.8) as 



= - / (1 - s)(w^22i[yk-i,yk-i^s{yk-yk-i)]{{yk-yk-i): {yk-yk-i),i^k-i) 
-^W^222[yk-i,yk-i^s{yk-yk-i)] 



{{yk-yk-i), (^/c-^fe-i), V^/c-i + s{iljk-i^k-i)) 



^2W^22[yk-i^yk-i^s{yk-yk-i)]{{yk-yk-i)^ (V^/c-V^/c-i))) ds 




for Vk = ^^^^ ^ and Wk = ^ . Now, taking into account the smoothness of the 

geodesic, the smoothness of W, and gy{v, w) = ^W,22[^, ^](^, which implies 



(5.10) 
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we finally achieve 
1 ^ 

= -^r(w^22i[yk-i,yk-i]{vk,Vk,'ilJk-i) + W^222[yk-i,yk-i]{vk,vk,'4^k-i) 

k—l \ 

-^2W^22[yk-i,yk-i]{vk,Wk)j +Err 

K 

= (2^2/^-1 K,^/c) + (^{Dygy^_J{iljk-i)Y^k^Vk)^ +Err (5.11) 

k=l 

for Err = 0(r2 5:j|^fc||3^||V;fc-i||v) + 0(r2 5:j|^fc||2^||^fc||v). From Theorem 5.6 we 
obtain the following estimates for the error, 

Err = 0(r2) (^^i llV'fc-i l|v + Ef=i Ikfcllv) , (5.12) 

Err = 0(r) (\/rEf=i ll^fc-i||^ + \/rj:ti Ik^llv) • (5-13) 

Based on these preliminaries we now combine the Euler-Lagrange equation (5.7) for the 
continuous geodesic path and (5.11) (derived from the Euler-Lagrange equation (5.8) for 
the discrete geodesic path) to obtain an equation for the discretization error. To this end, 
we consider a piecewise polygonal function ip^ : [0, 1] ^ V with i/jrikr) = G V for 
/c = 0, . . . , i^T and ipo = ipK = 0. Using this notation one easily verifies the identity 

= 2gy^t) {yr {t) -y{t), <P{t) - ijjr {t) ) dt 
-/o'2^,(,)(^(t),^,(t))+I),^,(,)(V^,(t))(^(t),^(t))dt 

2 [l^k-i)r 9yr it) {yr (t) , A (t) ) dt - Tgy^_^ {vf, , ^/c ) ) 

+2 Jq gy(t) {yr {t) , {t) ) - gy^ (t) {yr {t) , {t) ) dt 

+Ef=i (/(t-i)r Dygy^(t){i^k-i){yr{t),yr{t)) dt-TDygy^_^ (V^/c-l)K, ^/c)) 

(/(fc-l)r (^) ( V^r (^) ) (^/c , ^/c ) " ^jy^j/r (^) ( V^^-l ) {l/r {t) , (t))dt) 

+ /o^ {Dygyit){i^r{t))-Dygy^^t){i^r{t))) {yr{t) ,yr{t)) dt 

^IoDygy^t){^r{t)){Mt).y{t)-yr{t))^Dygy^t^{Mt)) 
=: I - II ^ III ^ IV ^VI ^ VII + VIII + /X (5.14) 

with v{t) = y{t), w{t) = ipr{t)- Now, we choose (}) = y^ — y and ijjj- = yr —^ry^ where X^- is 
the piecewise affine Lagrangian interpolation in time with Xrip{kr) = ip{kr). Let us mention 
here, that we will reuse (5.14) with different test functions later in the context of a pointwise 
error estimate in the proof of Theorem 6.1. Due to the uniform coercivity of the metric g the 
left-hand side can be estimated from below by 2c*\\y — ^r|li2((o i) V)- different terms 
on the right-hand side of (5.14) are estimated as follows: 

(I) By the uniform boundedness of the metric, | / | can be estimated by 

2C*\\y - ^r||L2((0,l);V)||^ " (^r^)' ||l2((0,1);V) • 

(II) Due to (5.7) the term // vanishes. 
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(III) Taking into account (5.11) and (5.13), | /// | is bounded by Cr\\ilJr\\w^''^{{o,i);V)- 

(IV) Due to the smoothness of the metric g and Theorem 5.6, the quadrature error in | IV \ 
is bounded by 0(r 2 ll^feHvll^r ||v) = C)(^||^r ||l2((o,i);V)) (compare (5.13)). 

(V) By an analogous argument and the boundedness of i/r in L^((0, 1); Y) the term | V \ 
can be estimated by ||l2((o,i);Y) ll^r ||l-((o,i);V) ||^r||L2((o,i);V), where we 

know from Theorem 5.6 that ||^r ||l°°((o,i);V) = ^(1)- 
(VI) Now using the smoothness of Dyg we can bound | VI \ by 

^^r'Ef=llKII^II«fe||Y||^fe-l||Y = O(r||^.|U.((0,l);Y)). 

(VII) Term | VII \ is bounded by Cr||?/'T-||/,2((o,i);Y) due to the uniform bound for ||^r||v- 
(VIII) Now using the smoothness of Dyg and the boundedness of ||^r||v we can bound 

I VIII I by C\\y - ||l2((0,1);Y) ll^r ||l2((0,1);Y). 

(IX) Finally, we obtain the estimate \IX\< ||l2((o,i);Y) ||^ - ||l2((o,i);V). 
Altogether, for some constant C we obtain the following estimate for the error, 

- ^r||i2((o,l);V) < - ^r||L2((0,l);V)||^ " (^r^)' ||l2((0,1);V) 

+Cr||V^^||iyi,2((o,i);V) 

-^C\\y - ^r||L2((0,l);Y)||^r||L2((0,l);V) 

+Cr||V^r||w^i'2((o,i);Y) 

^C\\y - ^r||L2((0,l);Y)||V^r||L2((0,l);Y) 
+C'||V^r||L2((0,l);Y)||^ - ^r||L2((0,l);V) • 

Finally, using Poincare's inequality \\y - ||l2((o,i);V) < 11^ - ||l2((o,i);V) as well as 
Young's inequality and the compact embedding V ^ Y, we derive in a straightforward way 

ll^-^r||i2((o,l);V) <C (t'^ ^ ll^-^r^||^i,2((o,l);V) + 11^ " ||i2((o,i);Y)) 

for some C > 0. Taking into account the classical interpolation estimate 

\\y - X^^||l2((0,1);V) + ^11^ - (^r^)1|L2((0,l);V) < Ct'^ 

as well as \\y — ||l2((o,i);Y) 0. this concludes the proof. □ 

Remark 5.12. Note that in the continuous setting, geodesies can also be obtained 
by minimizing the path length (2.1) for fixed end points ^(0), y{l). However, discrete min- 
imizers of the discrete path energy (2.1) are in general unrelated to continuous geodesies. 
As an example, think of a smooth curved manifold embedded in and let ^] be the 
squared Euclidean distance between y and y in the embedding space. Then given yo , yx ^ith 
yo + t{yi — Vo) for all t G (0, 1), minimizers ofL[-] for ?/o, yx fixed have the form 
{yo^yo^ • • • ^Vo^Vk^ • ' • ^yx) ' Obviously, the claim of Theorem 5.6 no longer holds, and the 
argument for Theorem 5.9 breaks down since minimizers of the discrete length L[-] are only 
bounded in 5V((0, l)]M) instead of W^^'^{{0, l);M). 

6. Convergence of discrete logarithm, exponential, and parallel transport. In this 
section we discuss the limit behavior of the discrete operators. At first we investigate the 
convergence of the discrete logarithm, which can be formulated as an derivative error 
estimate for discrete variational solutions of an elliptic problem in iy^'^((0, 1); A^). This 
will become apparent in the proof of the following theorem. 

Theorem 6.1 (Convergence of discrete logarithm). Let y^y ^ M. Under the hy- 
potheses (HI) and (H2), if the continuous and discrete geodesies between y^y are unique, 
K(^LOG) y log y weakly in V (and thus strongly inY) as K ^ oo. 
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Proof. As usual, abbreviate r = 1/K, and let ?/(t)^^[o,i] and ?/T-(^)tG[o,i] be the con- 
tinuous and the interpolated discrete geodesic between y and y. Let us denote by <l> G 
^0 '^((0, 1); Hom(V^ V)) the weak solution of 

r2^,(,)(i(t),i(t)t))dt= (t),^) (6.1) 

JO ^ /V',V 

for all z e 1^0 '^((0, 1); V), t) G V^ Then we obtain 

K(iLOG)^^ - log, y)^^^^ = ^-M) _ ^(0))^^^^ 

= /2^,(,)(^,(t)-^(t),^(t)t))dt + 0(r||t)||vO. (6.2) 

The solution of (6.1) for a given continuous geodesic path {y{t))te[o,i] can be computed 
explicitly. Indeed, we deduce from (6.1) the following two conditions, 

2^(%,)(^(t)t))) =0 on(0,r)U(r,l), (6.3) 
^%ir) (^(r + 0)t)-^(r-0)t)) =-^t) (6.4) 

for any t) G V^ where we use the notation : Y ^ V for the inverse Riesz isomorphism of 
the Hilbert space (V, gy). Integration of (6.3) leads to i>(t) = '^~^l^%(^o^^{0) on (0, r) and 
^{t) = on (r, 1). From this and (6.4) one obtains ^,(i)i>(l) - ^,(o)i>(0) = 

— Furthermore, due to the boundary conditions <l>(0) = <I>(1) = we achieve 

= l\t) dt ^j\t) dt = (yiy^il) dt) %io)m^(^iyy-^l^ d?j ^,(1)^(1X6.5) 

Now, with the notation G' = dt, G+ = jl dt, A = %(^o)^{0), and B = 

^,(i)i>(l), (6.4) and (6.5) form the linear system of equations 

G-A + G+B = 
which for 1 the identity has the unique solution 

A = -L(G-)-iG+(l + {G-)-'G+)-^ , B = -^{1 + (G-)-iG+)-i . 
zr zr 

Finally, we can evaluate ^{t) via integration of 

m = d^) ^y^il) A on (r, 1) . 
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Let US remark that in the trivial case of the constant metric gy{v, v) = (v, v)\- on a Hilbert 
space V (in which case geodesies are straight Hnes), <l> is a piecewise affine function with 
^{0) = <I>(1) = and ^{r) = for the Riesz isomorphism ^ : V ^ V. 

Next, we define as the Lagrangian interpolation X^^of^ with ^^{kr) = ^{kr) for 
k = 0, . . . , i^T. To show the claim, it remains to show the convergence of the right-hand side 
in (6.2). At first we proceed as in the proof of Theorem 5. 11. In (5.14) we choose ^ = 
and tpr = ^T^, which turns the left-hand side into the desired error representation. Hence, it 
remains to verify that the terms / to IX on the right-hand side of (5.14) vanish for r ^ 0. 
Indeed, using the same notation as in Theorem 5.1 1 we obtain the following estimates: 

(I) Due to the uniform boundedness of the metric, | / | can be bounded by 2C* \\y — 
^r||L2((o,i);V)||(*^ -^r*^*)'t5||L2((o,i);V) , and for t G {{k - l)r, /ct) one obtains 

< 'r I'^-Dr \Kil) - \\%it)^(i>\W = 0(r)||^(t>||v 

with t = for /c = 1 and i=l else. Using <i>(0) = 0{r-^) and <i>(l) = 0(1), we 
obtain ||(<l>-X^^)-t)||L2((o,i);V) = 0(r||t)||vO- 
(II) Due to (5.7) the term // vanishes as before. 

(III) By (5.1 1) and this time (5.12) the term | /// | can be estimated by 

C^||^r||L-((0,l);V) + Cr^||^l||v + Ct SUp \\Wk\W ^ 

k>l 

where ||i^;i||v < Cr-^||t)||v' andsup^^^ ll^feHv = 0(||t)||vO- 
Using ||V^r||L-((o,i);V) = 0(1), I /// I is bounded by Cr||t)||v'. 

(IV) Due to the smoothness of the metric g, the quadrature error in IV is bounded by 

0{t' J2 ll^fellYll^'fellvll^rllv) = 0{t^W,\W + TSUp \\Wk\W) = 0(t||d||v') • 

(V) Using Theorem 5.6, one obtains 

\V\<Cj^\\y-yA^\mW\\yr\W dt 

< C'l|y-yr||L~((0,l);Y)||^r||Li((0,l);V)l|yT||L°o((0,l);V) 
= O(||y-t/x||L~((0,l);Y)l|0||v')- 

(VI) Bythesmoothnessofi/weget|y/| < Cr^Ek IIV'fe-i||Y|kfc||^||vfe||Y = 0(T||t)||v') • 

(VII) Again due to the smoothness of g we achieve 

\VII\< Ct^wMv.WI, + Cj2T^wMv,r^ < Crllollv . 

k>l 

(VIII) For the term VIII one obtains the estimate 

I I < ^\\y - ^r||L2((0,l);Y)||V^r||L2((0,l);Y)||^r||ioo((o,l);V) 

= 0{\\y-yr\\LmO,l);Y)MW')- 

(IX) Finally, for the last term we get 

I /X I < C||?/^^||l2((0,1);Y)(||^||l-((0,1);V) + ||^r||L-((0,l);V))||^ " ||l2((0,1);V) 
= O(||^-^r||L2((0,l);V)||t)||v0- 
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Collecting all these estimates we obtain 

(v,K{^LOG)yy-\ogy y)^^^< C'(r+||7/-^,|U2((o,i);V) + lb-yr||L-((o,i);Y))||t)||v'. 

Using Theorems 5.10 and 5.1 1 the right hand side of this estimate converges to as i^T ^ oo 
for all t) G V^ which proves the claim. □ 

Next we consider the existence and convergence of the discrete exponential. This re- 
quires several preparations. 

Lemma 6.2 (Local uniqueness of LOG)). Under the hypotheses (HI) and (H2), there 
exists an e > such that for any yo^y2 ^ with y2 G B^{yo) = {y \ \\y — yo\\v ^ 
E[(yo7 '7^2)] • ^ is strictly convex with bounded coercive Hessian on B^{yQ), and 

(^L0G)^^(y2) is unique. 

Proof. For given e, consider ^0,^2 ^ with ||^2 — ^o||v < e. The Hessian of 
E[(^o, ^2)] at some?/i in the ball 5e(^o) is given by 2(W,22[^o, ^i]+W,ii[^i, ^2]) = 8^yo + 
0(e) with uniform constants, where we have used the smoothness of W and Lemma 5.7. For 
e small enough, this is bounded and coercive (independent of y^), and thus E[(?/o, ^2)] is 
strictly convex (cf. Theorem 5.8). This also implies uniqueness of (^LOG)^^ (7/2) (via The- 
orem 5.5 it is easy to see that the optimal yi cannot lie outside B^{yo)). □ 

Remark 6.3. The above result can be generalized to K -geodesies inside an e-ball. It is 
a refinement of Theorem 5.8 in the sense that it provides a uniform uniqueness radius e. Also, 
it implies the differentiability of {^LOG) ^^{y 2) with respect to y2 G B^{yo): The Hessian 

D'^E[{yo^ '7^2)] = 2(W^22[^07 '] + W 11 ['7^2]) is coercive so that by the implicit function 
theorem, 

D{^L0G)Jy2) = -(W,22[yo,^i] + W,ii[^i,y2])"'w,i2[^i,y2] 
= [^9yo + 0{e)]-\2gy, + 0(e)) = ^1 + 0(e) 

fory, :=yo^{lLOG)Jy2). 

Lemma 6.4. Fory2eBe{yo) the estimate ?/o+(|LOG)^^(?/2) = ^^+0(e^/2) holds. 

Proof Indeed, using >V[y,^] = dist^(?/, ^) + 0(dist^(y, ^)) = gyo{y-y,y-y)^0{e^) 
for y,y e B,{yo) and 2[W[^o, ^1(^2)] + W[^i(^2), ^2)]] = ^^0(^2 - ^0,^2 - ^0) + 0{e^) 
for yi{y2) = ^0 + (^L0G)^^(?/2) (by Theorem5.5), we get 

2^.0(^-1/1(^2), ^-^1(^2)) 

= 9yo {yo - yi {y2 ) , - (^2 ) ) + gyo (^2 - yi (^2 ) , ^2 - ^1 (^2 ) ) - ^ gyo (^2 -yo^y2-yo) 
= myo.yi{y2)]^myiiy2).y2] - (W[yo,^i(^2)]+w[yi(y2),^2]) + 0(e3) = 0{e^) . □ 

Lemma 6.5 (Local existence of (EXP^)). Under the hypotheses (HI) and (H2), there 
exists an e > such that (EXP^)^^('y) exists for any yo G M, v eV, with \\v\\\- < e. 

Proof. Choose e to be one third the value from Lemma 6.2. For yo ^ Ai, v ^ Y with 
ll^llv < e, define the operator Ty : Bs^iyo) M by Tyy2 = y2 ^ v - (^L0G)^^(?/2) • 
Without loss of generality we may assume e to be small enough such that by Lemma 6.4, 
Ty : Bseiyo) ^ Bseiyo)- Its derivative is given by 1 - I)(^LOG)^^ (^2) = ^1 + 0(e), 
hence Ty is a contraction, and by the Banach fixed point theorem has a unique fixed point 

{EXPX{v).n 

Remark 6.6. By the inverse function theorem, (EXP^)^^ is differentiable at v with 
\\v\W < e, and D{EXP\^{v) = [D{lL0G)J{EXP\{v))]-'=21 + 0{e). 

Lemma 6.7. For \\v\\v < e the estimate (EXP^) (w) = (yo + 2w) + 0{e^/'^) holds. 
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Proof. Let US abbreviate ?/i = ?/o+^ and 7/2(^1) = ^^^'^)yQ{^)- Then as in Lemma 6.4, 

\gy^{yQ^2v-y2{yi),yQ^2v-y2{yi)) 
= 9yo {yo -yi^yo-yi)^ 9yo (^2 (^1 )-yi,y2{yi)-yi)-l gyo {y2 {yi )-yo,y2{yi)-yo) 
= W[^o,^i] ^2(^1)]- (W[^o,^i] ^2(^1)]) + 0(6^) = 0(6^) .□ 

Remark 6.8. Using the above boundedness, our previous estimates can be further 
improved to \\{Wl^V\^{v) - (yo + 2t;)||v = 0{\\v\\\,) and ||(iL0G)^^(y2) - - 
yo)||v = 0(||t/2 - VoWlr) for \\v\\-y and ||t/2 - 2/o||v smaller than e. Indeed, if{yo,yi,y2) is 
a discrete geodesic with Vi = v = yi — yo and V2 = y2 — J/i. <^nd ifipi is a variation ofyi, 
then the optimality condition implies 

= W,2 [yo , yi] (^Ai ) + W,i [yi , ya] (^i ) 
= W,2[yi,yi](Vi)->V,2i[yi,yi](^/'i,yi -yo) + 0(||yi -yoll^ll^/'illv) 

+ >V,i[yi,yi](^/;i) + W,i2[yi,yi](V'i,y2-?/i) + 0(||y2 -yill^llV'illv) 
= 2y,,(^i, «i - V2) + 0((||«i + ||«2||^)||^i||v) , 

where we used Lemma 5.7. Using the coercivity of Qy^, the desired estimates can now easily 
be derived. 

Theorem 6.9 (Existence and convergence of (EXP^)). Let y : [0,1] M be a 
smooth geodesic. Under the hypotheses (HI) and (H2 ), (EXP^)^^q>j ( ) exists for K large 

enough, and for r = one obtains ||^(1) — (EXP^)^^q>j (^^) = 0(r) . 

Proof. To examine the convergence of the discrete exponential, we need to Hnearize 
the optimahty condition (5.8). Let (^o, • • • , ^k) be some discrete path in M and let Vk = 
{yk — yk-i)/^ for r = 1/K. As in the proof of Theorem 5. 11, for G V we find 

= r {1- ^) yv,222[yk-i,yk-i^tvk]{vk,Vk, ^i^k)^2yv^22[yk-i,yk-i^tvk]{vk, ^ipk) 

Jo L 



dt 







^ (1-r) - riSyki^m^A) - gyk-i{yk^A))^2Dygy^{tljk){Vk-i-l^Vk-H) 



J (2yV^222[yk-l,yk-l^rtVk]{Vk, ^Ipk.tVk) 

- 2 W,222 [yk , yk + rtVk-H ] {Vk-H , ^ V^/c , tVk-i-1 ) 

+ W,2212 [yk , ^/c , Vk-H ^ ^k , ) 

+ W,2222 [^fe, ^fe +^^^/c+l] (^/c+1, , V^/c, ^^/c+l)) dr 



+ W,222[^fe-l,^/c-l+^^/c](^fe,^fe, -V^/c) - W,222[^fe,y/c+^^/c+l](^/c+l,^/c+l, -^^/c) 
2 



dt 



r 



Dygy^iWiVk-i-l^Vk-H) 



T 

|2 



+ 0[(r(||^,||^+||^fe+i||^) + (r||^,||^+||^fe+i-^fc||v)(||^fc+i||v+|^ 
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assuming sufficient differentiabifity of W and using the simple Taylor expansion formula 

/[a2,62](c2,C2,6/)-/[ai,6i](ci,Ci,d) = /[ai,6i](c2-Ci,C2,c^)+/[ai,6i](ci,C2-Ci,d) 

+ /,i[a,6](c2,C2,c^, a2 - ai) + /,2[a, 6](c2, C2, c^, 62 - ^1) 

for a spatially dependent, differentiable trilinear form / and a = (1 — ^)ai + ^a2, b = 
(1 - 0^1 + ^^2 for some ^ G [0, 1]. 

Next, let us abbreviate tk = kr. The smooth geodesic y satisfies 

^9y{tk){'^k+i,i^k) - Qyitk-A'^k^i^k) . .. \ , ^/ II / II ^ 

= -2 +I)^^^(^^)(V^fc)K+i,^/c+i) + 0(r||V^fc||v) 

for = y^^^)-y^^^-^) and |0(r||V^,||v)| < Cr(||^'||v + ||l/||v(||^||v + r||i/||v))||V^.||v. 
Now assume yu+i •= {^^^^^^)y{()){y{^) / K) exists for some k < K (and thus also 
yoj ' ' ' J Vk)^ then the left-hand side of (6.6) is zero. Subtracting times (6.6) from the 
above equation yields 

O [(r + rll^fcll^ + r\\vk^i\\^ + {r\\vk\\v + ll^fe+i - ^fe l|v)(lbfc+i ||v + lbfe||v))||V^fe||v] 

T T 

- Dygy^{^l)k){vk+i,Vk+i) ^ Dygy^tk){'^k){uk+i,Uk+i) (6.7) 

Introducing Ck = yu — y{tk) as well as = ^ the first line on the right-hand side 

can be rewritten as 

r, 9yk{^k+l -Vk^i^k) - 9y{tk){^k+l -Uk.i^k) 
T 

^'2.Dygy^{vk){vk,^l)k) - 2Dygy^tk){'^k)(uk,^l^k) ^ 0{T{\\vk\\% + ||^/c||v)||^/c||v) 
=^9y\ ''^\'' .^k)^0[{\\e4^^^^ 

+ l|e^.llv(||^i.||v + \\vu\W) + r(||^,||^ + ||^.||^))||V^.||v] , 

and the second line is bounded by 0[(||efe||v||^fe+i||y + ||e^+i||v||^fe+i + ^/c+i||v)||V^/c||v]. 
Thus, exploiting the boundedness of y and its derivatives as well as 

Uk = y{tk) + rO{y) , Uk^i -Uk= ry{tk) + r^O(y) , 
\\vk\W < hk\W + \\el\W , H+1 - VkWv < hk+1 - UkWv + lle^+illv + ||e^||v , 

(6.7) can be rewritten as 

/ ek+i-Ck \ / \ 

V-^^J ~ UT+l|efe||v + (||e",||v + ||e",+il|v)(l + ||e^||v + K+illv)V ' 

where we have also used the uniform coerciveness of ^. If ||ej ||v < 1 for j = 0, . . . , /c + 1, 
then the above implies the existence of a constant C > such that with U = C (W), 
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which together with eo =0, ||ei ||v = 0(t) by classical arguments implies ||ej || v, ||v < 
nr, j = 0, . . . , /c + 1, for some k, > which is independent of k and K. 

Now choose ^ = K large enough such that Ck '= i^t'^ + rmax^^[o,i] ||l/(^)||v < 
e for e <C 1 from Lemma 6.5. Then by induction, (EXP^)(?/(0)/iir) = exists with 
Ikfcllv, ||e^||v < 1^^' Indeed, for /c = the situation is clear. Now if yj exists with 
l|ej||v,||ej||v < /^r for all j < A:, then - ^fe-2||v < (ll^/c-i ||v + ||e^_i ||v) r < 

(max^^[o,i] ||^(t)||v + /^r) r < < e so that y^ = (EXP^)^^_^(7//,_i - yk-2) exists by 

Lemma6.5 with yk - yu-i = Vk-i - yk-2 + 0(e^/^). Thus ||e/c||v, ||e^||v < 1, which via 
the above error estimates in turn implies ||e/e||v, ||e^||v < /^'r. □ 

Theorem 6.10 (Existence and convergence of P^^^...^2yo)- V • 1] ^ M be a 
smooth path, and let ( : [0, 1] Y be a parallel vector field along y. For G IN and 
r = let yk = y{kr), then under the hypotheses (HI) and (H2), the discrete parallel 

transport fulfills \\KPy^_y^{^) - C(l)||v = 0{r) . 

Proof. At first, we examine the residual resulting from the evaluation of the discrete 
counterpart of the parallel transport equation on the interpolated continuous vector field (. ( 
is defined to be parallel if ^( = holds for the covariant derivative ^ of ( along the curve 
X, where the covariant derivative is defined via 

9yit){{ft0{t).^) = 9y{t) (CW, ^) + 9y{t) iX{y(t). CW), ^) = (6.8) 

for all test vectors i/j e V and t e [0, 1]. Here the Christoffel tensor F : V x V ^ V is 
defined by gy{T{v, w), V^) = ^ [{Dygy){w){ilj, v) + {Dygy){v){ilj, w) - {Dygy){ilj){v, w)] 
for all v^w^i/j G V. Let us abbreviate tk = kr, (k = C(^fc)' = . Ap- 

plying gy{v^w) = ^yV^22[y,y]{v,w) (Lemma 5.7) and (5.10) as well as the approximations 

Citk-i) = ^^^1^ + 0(r||C||v), y{tk-i) = Vk^ 0(r||i/||v), (6.8) turns into 



O 



r{ sup ||CW||v+ sup \\y{t)\W sup ||CWIIv)||V^||v 

■ te[o,i] te[o,i] te[o,i] 



= W22[yk-uyk-i]{^-^'^^.^) 

+ ^ (yV,221 [yk-l,yk-l]{i^, Vk, Ck-l) + W,222[^fc-l, ^/c-l](V^, Vk, Ck-l) 

+W,22l[yfc-l,^/c-l](V^,C/e-l,^/e) + W,222 , ^/c-l] (V^, Cfc-l , ^/e) 

-W,22i [yk-i,yk-i] {vk, C/c-1, V^) - >V,222 [yk-i,yk-i] {vk, C/c-i, 

= W,22[yfc-l,yfe-l](^^^,V^) 

+ ^( - W,21l[^/c-l,^fe-l](V^,^fe,Cfe-l) - W,212[^fe-l,^fe-l](V^,^/c,C/c-l) 

+W,22i [yk-1, ^fe-i] (V^, , Vk) - >V,22i [yk-1, yk-i] {vk, C/c-i , V^)) 

= >V,22[^fc-l,^ife-l](^^^,^) 

(W,211 Vk^ Ck-l) + >V,22l[yfc-l, ^/c-l]K, C/c-1, V^)) 

r i/a-a-i /A W,ii2[yfe-l,^fe-l] + W,22l[y/c-l,yfe-l](Cfe-l,^/c,V^) 

= yy,22[yk-i,yk-i\[ — - — ,w) ^ 

where we have used Schwarz's theorem and the identity 

W,22i [y. y] + W,222 [y. y] = W,iii [y, y] + W,ii2 [y, y] 

= ->V,i2i [y, y] - W,i22 [y, ^] = - W,2ii - W,2i2 [y, ^] , 
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which can be derived by differentiating the identity W,22 ^] = VV,ii ^] = — W,i2 ^] = 
— W,2i [y-, y] from Lemma 5.7 with respect to y. Due to the smoothness of C and y we achieve 
the estimate 

- ^ (w,ii2 [yk-i.yk-i] + W,22i [yk-i.yk-i]) (C/c-i, ^/c, V^) = 0(r ||V^||v) (6.9) 

Next, we derive an estimate for the residual obtained when evaluating the discrete coun- 
terpart of the parallel transport equation on the discrete solution. Thus, let us abbreviate 
^/c = P2yfc,...,yo (C(0)/^) and denote the center of the /c* parallelogram in Schild's ladder by 
y^ . and y^ satisfy for all G V the two nonlinear equations 

y\^a[yk-i^ik-i.yim^y\^Aylykm = 0, (6.10) 
W^2[yk-i.yim^W^i[yl,yk^^km = 0, (6.11) 

which for the moment we assume to be uniquely solvable. Upon Taylor expansion about 
(?/£, y^) analogously to (5.9), one obtains 

y\^[yk-iHk-i.yl]= I {i-s)w^ii[yl^s{yk-iHk-i-yl).yl] 
Jo 

ivk-i^^k-i-yl^yk-i^^k-i-yl) ds, 



My"k.yk] = I (I - s)W^22[yl.yl^s{yk-yl)]{yk-yl.yk-yl) 
Jo 



ds. 



Differentiating these expressions, the first equation (6.10) defining the discrete parallel trans- 
port turns into 

\l - s) [ - 2W,n [yl, yl]{i>, yu-x +6-1 - yl) 

-2 [ W,iii[yl+rs{yk-i+^k-i-yi),yl]{^P,yk-i+^k-i-yl,s{yk-i+^k-i-y'i))dr 
Jo 

+>V,iii [y^+s(yfe_i +Cfe_i -y^), t/^](yfe_i +Cfe_i -y^, yk-i+^k-i-yt, (1 - sW 

+>v,ii2bfc+s(yfc-i +Cfc-i -yD. vlKyk-i +Cfc-i -yl, yk-i +6-1 ^/') 
-2-W,22[ylyl](^,yk-yl) 

-2 Wfi22[yl,yl+rs{yk - yl)]{^,yk - yl,s{yk - yl))dr 
Jo 

+yv,22i[yk,yk+s{yk - yk)]iyk - yk,yk - yti') 
+yv,222[yk,yk+4yk - vDMyk - yl^yk - yh{i - s)^/') ds 

-- -w,iibfe>2/fc](^,yfc-i+6-i -yfc) - w,22b£,y£](V',yfc -2/fe) 
+ \W,ix2[yl,yl]{yk-i+£.k-i - yl,yk-i+5.k-i - ylA) 
+ |>V,22i bfc, yfc] (yfe - t/fc, yfe - t/fc, ^/') 
+0[{\\yk-iHk-i-yl\\'^+\\yk-yl\\\ 



Here we used that for a spatially differentiable trilinear form / the relation —2/ [a, a 

rsv\{s'il),v,v) + /[a, a + ~ ^j'^^'^^'^) = ~ Ss)?/^, 'U, 'u) + 0(||'u||^ 
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holds for < r, 5 < 1, which after multipHcation with (1 — s) integrates up to 0(||v|p 
Likewise, the second equation (6.1 1) turns into 

= -yv,ii[yl, t/^](V', Vk-i - vl) - yv,22[yl y^KV', VkHk - vl) 
+|vy,ii2[yfe, yk\{yk-i -vh vk-i-Vk, ^) 
+\W,2-2i[yl,yl]{yk+£.k - yl^Vk+^k - ylA) 
+0 [{hk-i - yl\\\+\\ykHk - yl\\\mW] 

Subtracting the second from the first equation, dividing by r^, and using W n = W,22, one 
obtains 

H\}^MM.yk\{^.^-^='^^^^''^^^ 

+0 [{\\yk-i - yl\%+\\yk - y^ll^+IICfc-ill^+liail^)ll^llv] • (6.12) 

Finally, we derive an equation for the error propagation and subtract (6.9) from (6.12). Intro- 
ducing the error Ck = K^k — C/c = ^ — C/c. we arrive at 

WMvk, ytr-""^,^) = \ (>v,ii2bL yl] + w,22i [yl, yl]) {ek-i,Vk, ^) 



+ 



[t + \\yk-i-yU'v + hk-ytrv + Uk-ifv + liail^ + M-yk-iWvf-^h 



(6.13) 



+ U - yk-i\W\\Ck-i\W\H\W + ^||a-i||v||6-i - 2{yt - ^^^)||v 
+ 'rUk - ik-i\W\\vk\W + ^liailvlia - 2(y^ - ^^^^)||v] IIVIIv 

Applying the boundedness of Vk, Ck, C the uniform estimates 

11^ - S^llv < lia - a-illv + llefcllv + ||efc-i||v 
= 0(t + ||efe||v + ||efc-i||v) 
||^||v<||al|v + ||efc||v 
Ibfc - 2/fcllv = 0{\\yk-i + 6-1 - yfcllv) 

= TO(||«fe||/f + ||a-i||v + ||efe-i||v) 

hk - yk-i\W = 0{\\yk + ^k - yfe-i||v) 

= rO(||t;fc||v + ||ailv + ||efc||v) 
116-1 - 2{yl - y-^^)\W = 0{\\yk-^ + 6-1 - ykfv) 

= r^O{{\\vk\W + \\Ck-i\W + \\ek-i\W?) 
Uk - 2{yl - ^^^±f^)||v = 0{\\yk-, - yk - ik\\lr) 

= r^O{{\\vk\\K + \\Ck\W + \\ek\W?) 

(the last four follow from Remark 6.8, assuming ||^/e_i||v, ||^^/c||v to be sufficiently small), 
we find 

2^^^ = 0[||efc-i||v + t(1 + \\ek-i \Wf + ||efc||v + t(1 + ||efc||v)'] ■ (6.14) 



As in the proof of Theorem 6.9, this implies the existence of a constant C > (independent 
of K) such that if ||e^||v <lforj<k<K, then ||efe||v < Cr. 
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We now inductively show the existence of ^/e with ||e/c||v ^ Cr. Choose K large 
enough such that C + sup;.^i^___^^ ||C/e-i||v + lb/e||v + ||C/e||v < ei^ for a given e <C 1 from 
Lemma 6.5. Clearly, exists, and eo = 0. Now assume the existence of ^j, j < k, with 
||ej||v ^ Cr. In particular, this implies ||^/e_i||v + T||v/e||v < e and thus by Lemmata 6.2 
and 6.5 the existence of and ^/e. Furthermore, Remark 6.8 ensures ||e/c||v ^ 1 such that 
the above estimates imply ||e/c||v ^ Cr. □ 

The parallel transport also converges if the interpolated smooth path (?/o, . . . , ^k) is 
replaced by another approximating path, for instance a discrete geodesic. 

Corollary 6. n . Let y : [0,1] ^ M be a smooth path and C : [0, 1] ^ V a parallel 
vector field along y. For K G ¥i andr = j^, let yk satisfy \\yk—y{]^^)\W < e, /c = 0, . . . , 
then under the hypotheses (HI) and (H2), WKVyj^^^^^^y^^^^) — C(l)||v = 0{r + e) . 

Proof. Following the previous proof we arrive at a slightly altered version of (6.13), 



-^M- ^(ti.-i)llv||a-i||v||^.||v + ^lia-illvlia-i - 2{yl - ^^^)||v 

+ ^lia - a-i||v||^.||v + ^liailvlia - 2{yl - ^^^)||v] ll^llv , 
where t/e = kr^n&Vk = y^^^^-y^^^-^\ Using -^(tfc-i)llv < - ^/c-i ||v + e as well 



as the estimates from the previous proof we obtain (6.14) with 0(e) added to the right-hand 
side, from which the claim follows by the same arguments. □ 

Corollary 6.12. Let ^ g V and j]-. yA^^^^bea smooth vector field. Define 
= iviv)^ v{y ~^ '^^))' ^^^^ under the hypotheses (HI) and (H2)we have || Vr6>(T7^^) — 
Vev\W = 0{r). 

Proof. Consider the continuous geodesic y{t) with y{0) = y and ^(0) = 0. We have 
\\y{r) - (y + re)\W = 0{r^)- Now define ^ = rr]{y + rO) and = ^ylroji ^s well 
as (i = r]{y{r)) and Co as the vector parallel transported from y{r) to y{0). Furthermore 
introduce the error Ck = — Ck, k = 0,1. Since P~+^5) y is defined via the same discrete 
Euler-Lagrange equations ((2.6) and (2.7) for k = 0, yo = y, and yi = y -\- rO) as Vyj^^e^y^ 
(6. 14) holds also here, only withO(||?/(r) — (?/+r6>)||v) = O(r^) added to the right-hand side 
(exactly as in the previous corollary). Using ei = O(r^), we immediately see eo = O(r^). 
The relation V eV = ^^^^ + 0(r) together with Co + 0{r'^) = ^^o = r]{y) + ^Vre{rr]^) 
now implies the claim. □ 

7. Application to shape spaces. In this section we will discuss the scope of applications 
of the presented convergence theory. In general, it can be utilized in a straightforward fashion 
for finite-dimensional manifolds and for those infinite-dimensional manifolds where motion 
paths are compact in state space and the metric can be rephrased properly on Lagrangian 
velocity fields. In what follows we exemplify this in the case of the discrete geodesic calculus 
on embedded manifolds from Section 4, on the parametrized finite-dimensional shape space 
of viscous-fluid objects from Section 3.1, and on the simplified shape space of viscous one- 
dimensional shells from Section 3.2. 

In its current form the convergence theory is not yet applicable to the flow of diffeomor- 
phism approach in Section 3.3. The difficulty lies in identifying the manifold M locally with 
an affine space + V of equivalent topology, in particular if shape variations are described 
in infinite dimensions by Eulerian instead of Lagrangian velocity fields. 



Cfc~Cfc-i 



r 



Ilv 



7.1. On embedded finite-dimensional manifolds (cf. Section 4). The most basic ex- 
ample has been given in Section 4 by a smooth, complete, (m — 1) -dimensional manifold M, 
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embedded in W^. Due to the assumptions on A^, a continuous geodesic exists between any 
two points Va^Vb ^ M.- Furthermore, we assume the embedding to be such that (5.5) holds 
for the function W which represents the squared extrinsic distance between two points on M . 

Even though in this example M. is not an affine space, the existence proof for discrete 
geodesies (Theorem 5.4) still applies if in it we replace Y hy M. (note that due to the fi- 
nite dimensions, the proof no longer has to rely on a compact embedding of M into some 
larger space). All subsequent results either do not make use of the manifold structure or 
are essentially local in the sense that their proofs only consider a neighborhood of a con- 
tinuous geodesic. This characteristic enables us to transfer the results onto the example 
at hand: Via a chart we locally identify the manifold M. with an open subset of IR^~^ 
and set Y = V = IR^~^. In more detail, all convergence proofs involve a continuous 
geodesic or smooth path (the limit object of convergence), around which a small compact 
tubular neighborhood U can be identified with a compact subset ^(/7) of IR^~^ (where 
(j) : U ^ IR^~^ is the parametrization). The induced metric on 0(/7) is clearly equiva- 
lent to the Euclidean metric. Hence, if we consider (t)(U) C Y instead of A1, we are in the 
setting of Sections 5 and 6, and all convergence results apply. Now, convergence of paths yr 
inl^^'^((0, 1) ;(/)([/)) and of vectors ^Vj- G V against vectors v in a tangent space T^^^) (/)([/) 
immediately implies the same rate convergence of geodesies (j)~^{y^) in iy^'^((0, 1); U) and 
of ^{(})~^{(j){y) + Vr) — y) ^ against D(})~^v G TyU so that discrete geodesies, loga- 
rithm, exponential map, and parallel transport converge as well in this sense. From the proofs 
it should be clear that the restriction to U does not impair the results. 

7.2. On the parametrized finite-dimensional shape space of viscous-fiuid objects 
(c/. Section 3.1). If we consider viscous fluid objects being represented by spline-parame- 
trized boundary contours, each shape can be represented (uniquely up to rigid body motions) 
by an n-tuple (^i, . . . , ^n) of control points in R'^ (d = 2 or d = 3). First we quotient 
out the rigid body motion invariance. To this end we fix a reference configuration yA = 
(gi, . . . , Qn) G MJ^^ with Ylk^i Qk = ^ and only consider those shapes y = (g'l, . . . , g^) that 
can be reached by a displacement of yA with zero average linear and angular momentum, i.e. 
y eM = yA^yforV = {{qi,...,qn) eR^^ \ ELi = and ELi Qk ^ Qk = 0} . 
Also, we define Y = V. Then Y is a /c-dimensional subspace of (k = 2{n — 1) — 1 
for (i = 2 and k = 3{n — 2) for d = 3), of which the admissible shapes actually only make 
up a closed subset, i.e. a manifold with boundary, since shapes with self-intersection of the 
shape contour are not allowed. In fact, we shall restrict the admissible shapes to those whose 
boundaries satisfy a uniform cone condition, i.e. each point on the shape boundary must be 
the tip of two open cones with fixed opening angle a and height r, one contained inside and 
one outside the shape. In this context, the metric is given by Qyiy^ y) = Jo^y^ + 
/itr(e[i;]^) dx for = (g'l, . . . , Qn) G Y, = (gi, . . . , Qn) ^ V, where 0[y] is the interior 
of the spline described by y, and v : 0[y] ^ R^ is the velocity field minimizing the right 
hand side with Dirichlet boundary data prescribed by the variation of y and thus of all points 
on dO[y]. Using Korn's inequality, it is a classical result in the direct method of the calculus 
of variations that a minimizing v exists. This v depends linearly on y, hence gy is indeed 
a bounded quadratic functional. Again by Korn's inequality, gy{y^y) = implies v{x) = 
vq + Sx for a constant velocity vq and a skew- symmetric matrix S. However, vq ^ or 
S contradict y = {v{qi)^ . . . ,v{qn)) G V, which implies that gy is positive definite, 
and the shape space indeed carries the structure of a Riemannian manifold. Furthermore, 
by the existence and regularity theory of linearized elasticity [ ], using the implicit function 
theorem, the metric gy and the functional W actually smoothly depend on y with gy = 
\y^,22[y-,y\ so that (2.4) holds by Lemma 5.7. Furthermore, by the direct method in the 
calculus of variations we see that for points yA and ys which are sufficiently close there 
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exists an energy minimizing path (the imposed cone condition for the shapes renders the 
existence of smooth geodesic paths for arbitrary end shapes yA and ys unclear). Now, from 
the coercivity property (5.5), we know that any discrete iiT-path between two points yAjl/B ^ 
M with smaller energy than E[(^^, . . . , ^a, ^b)] has to lie completely within a radius of 
K'y~^{E[{yA, . . . , yA^ yB)]/K) around yA- Hence, in the proof of Theorem 5.4 we may 
restrict to that compact region. This yields existence of a discrete iC-geodesic for any K. 
Theorem 5.6 now implies that all K-geodesics between yA and ys stay inside the same closed 
bounded and thus compact region. Hence, during all our convergence analysis we may restrict 
to a compact subset of = V = Y on which the Euclidean metric and g are equivalent, 
hence all convergence results hold. 

7.3. On the shape space of viscous shells (cf. Section 3.2). For the simplified model 
of one-dimensional shells (3.3), (3.4) let us consider for < a < ^ the space Y = {?/ G 
C^'"(IR;IR^) \ y{s + 1) = y{s)} of C^'"-smooth closed curves represented by periodic 
parametrizations. Furthermore, let V = G W^'^(IR; R^) | v{s + 1) = v{s)^ v{s) ds = 
^ v{s) AyA{s)ds = 0} for some fixed G Y be the space of velocity fields with zero 
average speed and angular momentum. Furthermore, we assume that \yA,s\ ^ c > 0. Then, 
we identify the manifold Ai of one-dimensional shells locally with an open subset of + V. 
By Sobolev embedding results V is compactly embedded in Y, and V is a reflexive, sepa- 
rable space. Given v G i^^((0, 1); V) with \\v{t^ •)llvv2'2((o i) dt < e for sufficiently 

small e we obtain from y{t^ s) = yA{s) + Jq v(r, s) dr that | < |y,s(^, s)\ < C for the path 
y G i^^((0, l);M) with v = y. Furthermore, for gy defined in (3.4) gy{v^v) = implies 
^' = 0. From this we deduce that in a sufficiently small neighborhood of yA the metric g is 
uniformly bounded and uniformly coercive on V in the sense of (5. 1) and that y ^ gy and W 
defined in (3.3) are smooth. Furthermore, the coercivity assumption (5.5) for W is fulfilled. 
Hence, for yB — yA sufficiently small the direct method in the calculus of variations used 
in the proofs of Theorem 5.3 and Theorem 5.4 can be applied to establish the existence of a 
continuous and a discrete geodesic. Furthermore, Theorem 5.6 implies that a iC-geodesic be- 
tween yA and ys stays inside the same bounded region in + V. As above, the convergence 
analysis can be restricted to such a small neighborhood of yA- 

For the discrete approach on triangular meshes of fixed topology we describe a shape via 
the vertex positions and choose Y = V = R^'^, where n is the number of nodes. Due to the 
smoothness of the energy W (3.2) our convergence theory directly applies even though the 
geometric interpretation of the associated metric gy = ^ W,22 [y^ v] is still open. 

Acknowledgements. The authors thank Behrend Heeren for providing Figure 3.4. 
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